Moving from Mayavi to Matplotlib

A long standing dependency of MotionClouds is MayaVi. While powerful, it is tedious to compile and may discourage new users. We are trying here to show some attempts to do the same with matplotlib.

In [2]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Testing 3D plots in pylab

In [3]:
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
cset = ax.contour(X, Y, Z, cmap=cm.coolwarm)
ax.clabel(cset, fontsize=9, inline=1)
  Hide output
In [4]:
X, Y, Z
  Hide output
Out[4]:
(array([[-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        ..., 
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5],
        [-30. , -29.5, -29. , ...,  28.5,  29. ,  29.5]]),
 array([[-30. , -30. , -30. , ..., -30. , -30. , -30. ],
        [-29.5, -29.5, -29.5, ..., -29.5, -29.5, -29.5],
        [-29. , -29. , -29. , ..., -29. , -29. , -29. ],
        ..., 
        [ 28.5,  28.5,  28.5, ...,  28.5,  28.5,  28.5],
        [ 29. ,  29. ,  29. , ...,  29. ,  29. ,  29. ],
        [ 29.5,  29.5,  29.5, ...,  29.5,  29.5,  29.5]]),
 array([[-0.01 , -0.011, -0.013, ..., -0.015, -0.013, -0.011],
        [-0.011, -0.013, -0.015, ..., -0.018, -0.015, -0.013],
        [-0.013, -0.015, -0.018, ..., -0.02 , -0.018, -0.015],
        ..., 
        [-0.012, -0.014, -0.017, ...,  0.029,  0.03 ,  0.031],
        [-0.011, -0.013, -0.015, ...,  0.016,  0.017,  0.018],
        [-0.01 , -0.012, -0.014, ...,  0.007,  0.008,  0.009]]))
In [6]:
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np

n_angles = 36
n_radii = 8

# An array of radii
# Does not include radius r=0, this is to eliminate duplicate points
radii = np.linspace(0.125, 1.0, n_radii)

# An array of angles
angles = np.linspace(0, 2*np.pi, n_angles, endpoint=False)

# Repeat all angles for each radius
angles = np.repeat(angles[...,np.newaxis], n_radii, axis=1)

# Convert polar (radii, angles) coords to cartesian (x, y) coords
# (0, 0) is added here. There are no duplicate points in the (x, y) plane
x = np.append(0, (radii*np.cos(angles)).flatten())
y = np.append(0, (radii*np.sin(angles)).flatten())

# Pringle surface
z = np.sin(-x*y)

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_trisurf(x, y, z, cmap=cm.jet, linewidth=0.2, alpha=.6)
  Hide output
Out[6]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x10a16d590>
In [10]:
"""
Contour plots of unstructured triangular grids.
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as tri
import numpy as np
import math

# First create the x and y coordinates of the points.
n_angles = 48
n_radii = 8
min_radius = 0.25
radii = np.linspace(min_radius, 0.95, n_radii)

angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
angles = np.repeat(angles[...,np.newaxis], n_radii, axis=1)
angles[:,1::2] += math.pi/n_angles

x = (radii*np.cos(angles)).flatten()
y = (radii*np.sin(angles)).flatten()
z = (np.cos(radii)*np.cos(angles*3.0)).flatten()

# Create a custom triangulation
triang = tri.Triangulation(x, y)

# Mask off unwanted triangles.
#xmid = x[triang.triangles].mean(axis=1)
#ymid = y[triang.triangles].mean(axis=1)
#mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
#triang.set_mask(mask)

plt.figure()
plt.gca(projection='3d')
plt.tricontour(triang, z)
  Hide output
Out[10]:
<matplotlib.tri.tricontour.TriContourSet instance at 0x10bdd23f8>
In [12]:
from __future__ import print_function
"""
A very simple 'animation' of a 3D plot
"""
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
import time

def generate(X, Y, phi):
    R = 1 - np.sqrt(X**2 + Y**2)
    return np.cos(2 * np.pi * X + phi) * R

plt.ion()
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

xs = np.linspace(-1, 1, 50)
ys = np.linspace(-1, 1, 50)
X, Y = np.meshgrid(xs, ys)
Z = generate(X, Y, 0.0)

wframe = None
tstart = time.time()
for phi in np.linspace(0, 360 / 2 / np.pi, 100):

    oldcol = wframe

    Z = generate(X, Y, phi)
    wframe = ax.plot_wireframe(X, Y, Z, rstride=2, cstride=2)

    # Remove old line collection before drawing
    if oldcol is not None:
        ax.collections.remove(oldcol)

    plt.pause(.001)

print ('FPS: %f' % (100 / (time.time() - tstart)))
  Hide output
FPS: 175.698996

However, there is nothing really easy yet to produce a contour plot of a scalar field, as mayavi does very well.

Testing the utility functions of Motion Clouds

Motion Clouds utilities

Here, we test some of the utilities that are delivered with the MotionClouds package.

In [1]:
import os, sys
sys.path.append('..')
import MotionClouds as mc
  Hide output

progressbar

Wrapper around pyprind.

In [2]:
### PERCENTAGE INDICATOR EXAMPLE ###
n = 10000
my_prperc = mc.progressbar.ProgPercent(n) # 1) init. with number of iterations
for i in range(n):
    # do some computation
    my_prperc.update() 
  Hide output
[100 %] elapsed[sec]: 8.578
Total time elapsed: 8.579 sec

In [9]:
name = 'color'
filename = os.path.join(mc.figpath, name)
print filename, mc.check_if_anim_exist(filename), mc.recompute
%bookmark root
%cd -b root
%cd..
  Hide output
results/color True False
(bookmark:root) -> /Users/lolo/Dropbox/science/MotionClouds/posts
/Users/lolo/Dropbox/science/MotionClouds/posts
/Users/lolo/Dropbox/science/MotionClouds

In [10]:
print mc.MAYAVI
mc.MAYAVI = 'Import'
mc.import_mayavi()
print mc.MAYAVI
  Hide output
Could not import Mayavi
Could not import Mayavi
Could not import Mayavi

exporting to various formats

In [13]:
name = 'export'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.rectif(mc.random_cloud(mc.envelope_gabor(fx, fy, ft)))
for vext in mc.SUPPORTED_FORMATS:
    mc.anim_save(z, os.path.join(mc.figpath, name), display=False, vext=vext)
  Hide output
[6350 %] elapsed[sec]: 0.903 | ETA[sec]: -0.888 
Saving sequence results/export.mpg
Title: 
                      Started: 06/18/2014 17:10:26
                      Finished: 06/18/2014 17:10:26
                      Total time elapsed: 0.000 sec
[  0 %] elapsed[sec]: 0.000

Saving sequence results/export.mp4
[6350 %] elapsed[sec]: 0.836 | ETA[sec]: -0.823 

Title: 
                      Started: 06/18/2014 17:10:27
                      Finished: 06/18/2014 17:10:27
                      Total time elapsed: 0.000 sec
[  0 %] elapsed[sec]: 0.000

Saving sequence results/export.gif
[6350 %] elapsed[sec]: 0.871 | ETA[sec]: -0.857 

Title: 
                      Started: 06/18/2014 17:10:28
                      Finished: 06/18/2014 17:10:28
                      Total time elapsed: 0.000 sec
[  0 %] elapsed[sec]: 0.000

Saving sequence results/export.zip
[6350 %] elapsed[sec]: 0.765 | ETA[sec]: -0.753 

Title: 
                      Started: 06/18/2014 17:10:29
                      Finished: 06/18/2014 17:10:30
                      Total time elapsed: 0.000 sec
[  0 %] elapsed[sec]: 0.000

Saving sequence results/export.mkv
[6350 %] elapsed[sec]: 0.808 | ETA[sec]: -0.796 

Title: 
                      Started: 06/18/2014 17:10:34
                      Finished: 06/18/2014 17:10:35
                      Total time elapsed: 0.000 sec

showing a video

In [12]:
mc.in_show_video('color')
  Hide output
color

In []:
 
  Hide output
In []:
 
  Hide output

Testing the different components of the enveloppe

Motion Clouds: testing components of the envelope

In [2]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output
In [3]:
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
#mc.recompute = True
mc.notebook = True
  Hide output

Testing the speed

Here the link to the test page for the component Speed

In [5]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
name = 'speed'
z = color*mc.envelope_speed(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
speed

Exploring the orientation component of the envelope around a grating.

Here the link to the test page for the component Grating

In [6]:
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
color = mc.envelope_color(fx, fy, ft)
name = 'grating'
z = color * mc.envelope_gabor(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
grating

Here the link to the test page for the component Radial

In [8]:
B_theta = np.pi/8.
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
name = 'radial'
mc.figures_MC(fx, fy, ft, name, B_theta=B_theta)
verbose = False
mc.show_video2(name)
  Hide output
radial

Testing the color

Here the link to the test page for the component Color

In [10]:
name = 'color'
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
z = mc.envelope_color(fx, fy, ft)
mc.figures(z, name)
mc.show_video2(name)
  Hide output
color

Testing competing Motion Clouds with psychopy

Analysis of a Discrimination Experiment

In the psychopy_competition.py script, we implemented an experiment to test whether one could discriminate the dominating motion in a sum of two motion clouds in opposite directions.

Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.

In [27]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Analysis, version 1

In a first version of the experiment, we only stored the results in a numpy file:

In [28]:
!ls ../data/competing_v1_*npy
  Hide output
../data/competing_v1_bruno_Dec_14_1210.npy  ../data/competing_v1_lup_Dec_12_1013.npy  ../data/competing_v1_meduz_Dec_14_1204.npy
../data/competing_v1_lup_Dec_12_1003.npy    ../data/competing_v1_lup_Dec_14_1201.npy

In [29]:
import glob
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    print 'experiment ', fn, ', # trials=', results.shape[1]
    ax.scatter(results[1, :], results[0, :])
_ = ax.axis([0., 1., -1.1, 1.1])
_ = ax.set_xlabel('contrast')
_ = ax.set_ylabel('perceived direction')
  Hide output
experiment  ../data/competing_v1_bruno_Dec_14_1210.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_12_1003.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_12_1013.npy , # trials= 50
experiment  ../data/competing_v1_lup_Dec_14_1201.npy , # trials= 10
experiment  ../data/competing_v1_meduz_Dec_14_1204.npy , # trials= 50

In [30]:
alpha = .3
fig = plt.figure(figsize=(12,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
data = []
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    data_ = np.empty(results.shape)
    # lower stronger, detected lower = CORRECT is 1
    ax1.scatter(results[1, results[1, :]<.5], 
                1.*(results[0, results[1, :]<.5]==-1), 
                c='r', alpha=alpha)
    # copying data: contrast (from .5 to 1), correct 
    data_[0, results[1, :]<.5] = 1. - results[1, results[1, :]<.5]
    data_[1, results[1, :]<.5] = 1.*(results[0, results[1, :]<.5]==-1)
    # upper stronger, detected lower = INCORRECT is 1
    ax1.scatter(results[1, results[1, :]>.5], 
                1.*(results[0, results[1, :]>.5]==-1), 
                c='b', alpha=alpha)
    # lower stronger, detected upper = INCORRECT is 1
    ax2.scatter(results[1, results[1, :]<.5], 
                1.*(results[0, results[1, :]<.5]==1), 
                c='b', alpha=alpha, marker='x')
    # upper stronger, detected upper = CORRECT is 1
    ax2.scatter(results[1, results[1, :]>.5], 
                1.*(results[0, results[1, :]>.5]==1), 
                c='r', alpha=alpha, marker='x')
    # copying data: contrast (from .5 to 1), correct 
    data_[0, results[1, :]>=.5] = results[1, results[1, :]>=.5]
    data_[1, results[1, :]>=.5] = 1.*(results[0, results[1, :]>=.5]==1)
    data.append(data_)

for ax in [ax1, ax2]:
    ax.axis([0., 1., -.1, 1.1])
    _ = ax.set_xlabel('contrast')
_ = ax1.set_ylabel('detected lower')
_ = ax1.legend(['lower', 'upper'], loc='right')
_ = ax2.set_ylabel('detected upper')
_ = ax2.legend(['lower', 'upper'], loc='right')
  Hide output

(note the subplots are equivalent up to a flip)

One could not fit psychometric curves for this 2AFC. See for instance https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3110332/

let's explore another way:

In [31]:
import numpy as np
import pylab
from scipy.optimize import curve_fit
 
def sigmoid(c, c0, k):
    y = 1 / (1 + np.exp(-k*(c-c0)))
    return y
 
  Hide output
In [33]:
for fn in glob.glob('../data/competing_v1_*npy'):
    results = np.load(fn)
    cdata, ydata = results[1, :], .5*results[0, : ]+.5
    pylab.plot(cdata, ydata, 'o')

    popt, pcov = curve_fit(sigmoid, cdata, ydata)
    print popt
    c = np.linspace(0, 1, 50)
    y = sigmoid(c, *popt)
    pylab.plot(c, y)
    
pylab.ylim(-.05, 1.05)
pylab.legend(loc='best')
  Hide output
[ 0.291  6.067]
[  0.345  13.492]
[  0.287  39.2  ]
[   0.637  160.796]
[  0.362  13.633]

Analysis, version 2

In the second version, we also stored some information about the experiment

In [41]:
from psychopy import misc
mean, slope = [], []
for fn in glob.glob('../data/competing_v2_*npy'):
    results = np.load(fn)
    data = misc.fromFile(fn.replace('npy', 'pickle'))
    print data
    cdata, ydata = results[1, :], .5*results[0, : ]+.5
    pylab.plot(cdata, ydata, 'o')

    popt, pcov = curve_fit(sigmoid, cdata, ydata)
    mean.append(popt[0])
    slope.append(popt[1])
    c = np.linspace(0, 1, 50)
    y = sigmoid(c, *popt)
    pylab.plot(c, y)

pylab.text(0.05, 0.8, 'mean : %0.2f , slope : %0.2f ' %(np.mean(mean), np.mean(slope)))
pylab.ylim(-.05, 1.05)
pylab.legend(loc='best')
  Hide output
{'N_X': 128, 'N_Y': 128, 'observer': u'jean', 'timeStr': 'Sep_03_1536', 'nTrials': 50, 'N_frame': 128, 'N_frame_total': 128, 'screen_width': 2560, 'screen_height': 1440}
{'N_X': 128, 'N_Y': 128, 'observer': u'laurent', 'timeStr': 'Sep_17_1522', 'nTrials': 50, 'N_frame': 128, 'N_frame_total': 128, 'screen_width': 2560, 'screen_height': 1440}

Testing discrimination between Motion Clouds with psychopy

Analysis of a discrimination experiment

In the psychopy_discrimination.py script, we implemented an experiment to test whether the overall shape of motion clouds could influence discrimination of speed.

Herein, we analyse the data that was collected over different sessions and try to draw some conclusions.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

The experiment

In [2]:
#%pycat psychopy_discrimination.py
  Hide output

Analysis - version 1

In a first version of the experiment, we only stored the results in a numpy file.

In [3]:
!ls  ../data/discriminating_*
  Hide output
../data/discriminating_john_Jan_23_1515.npy  ../data/discriminating_lup_Jan_24_0914.npy  ../data/discriminating_lup_Jan_24_0937.npy
../data/discriminating_lup_Jan_23_1401.npy   ../data/discriminating_lup_Jan_24_0919.npy  ../data/discriminating_v2_lup_Feb_07_1409.pickle
../data/discriminating_lup_Jan_23_1735.npy   ../data/discriminating_lup_Jan_24_0931.npy  ../data/discriminating_v2_lup_Feb_07_1434.pickle

In [4]:
!ls  ../data/
  Hide output
competing.pickle		    competing_v1_meduz_Dec_14_1204.npy	       discriminating_john_Jan_23_1515.npy  discriminating_lup_Jan_24_0931.npy
competing_v1_bruno_Dec_14_1210.npy  competing_v2_anonymous_Feb_07_1558.pickle  discriminating_lup_Jan_23_1401.npy   discriminating_lup_Jan_24_0937.npy
competing_v1_lup_Dec_12_1003.npy    competing_v2_anonymous_Sep_16_1746.pickle  discriminating_lup_Jan_23_1735.npy   discriminating_v2_lup_Feb_07_1409.pickle
competing_v1_lup_Dec_12_1013.npy    competing_v2_jean_Sep_03_1536.npy	       discriminating_lup_Jan_24_0914.npy   discriminating_v2_lup_Feb_07_1434.pickle
competing_v1_lup_Dec_14_1201.npy    competing_v2_jean_Sep_03_1536.pickle       discriminating_lup_Jan_24_0919.npy

In [5]:
import glob
for fn in glob.glob('../data/discriminating_*npy'):
    results = np.load(fn)
    print fn, results.shape
    print('analyzing results')
    print '# of trials :', np.abs(results).sum(axis=1)
    print 'average results: ', (results.sum(axis=1)/np.abs(results).sum(axis=1)*.5+.5)
  Hide output
../data/discriminating_john_Jan_23_1515.npy (2, 50)
analyzing results
# of trials : [ 50.     24.508]
average results:  [ 0.48  1.  ]
../data/discriminating_lup_Jan_23_1401.npy (2, 50)
analyzing results
# of trials : [ 50.     28.126]
average results:  [ 0.66  1.  ]
../data/discriminating_lup_Jan_23_1735.npy (3, 50)
analyzing results
# of trials : [  9.  14.  13.]
average results:  [ 1.  1.  1.]
../data/discriminating_lup_Jan_24_0914.npy (3, 50)
analyzing results
# of trials : [ 17.  21.  12.]
average results:  [ 0.647  0.857  1.   ]
../data/discriminating_lup_Jan_24_0919.npy (3, 50)
analyzing results
# of trials : [ 10.  25.  15.]
average results:  [ 0.7    0.32   0.533]
../data/discriminating_lup_Jan_24_0931.npy (7, 50)
analyzing results
# of trials : [  3.   4.   8.   8.   7.  12.   8.]
average results:  [ 0.667  1.     0.625  0.375  1.     0.167  0.125]
../data/discriminating_lup_Jan_24_0937.npy (7, 50)
analyzing results
# of trials : [  7.   5.   6.  12.  10.   2.   8.]
average results:  [ 0.857  0.8    0.833  0.417  0.2    1.     0.375]

In [6]:
%whos
  Hide output
Variable   Type       Data/Info
-------------------------------
fn         str        ../data/discriminating_lup_Jan_24_0937.npy
glob       module     <module 'glob' from '/usr<...>/lib/python2.7/glob.pyc'>
np         module     <module 'numpy' from '/us<...>ages/numpy/__init__.pyc'>
plt        module     <module 'matplotlib.pyplo<...>s/matplotlib/pyplot.pyc'>
pylab      module     <module 'pylab' from '/us<...>site-packages/pylab.pyc'>
results    ndarray    7x50: 350 elems, type `float64`, 2800 bytes

Another solution is to use:

In [7]:
#data= 1
#np.savez('toto', data=data, results=results)
#print np.load('toto.npz')['data']
  Hide output

or directly psychopy.misc:

In [8]:
from psychopy import misc
  Hide output
In [9]:
#info = misc.fromFile('../data/discriminating.pickle')
  Hide output
In [10]:
alphas = [-1., -.5, 0., 0.5, 1., 1.5, 2.]
fileName = '../data/discriminating_' + info['observer'] + '.pickle'
info['alphas'] = alphas
info['results'] = results
#numpy.savez(fileName, results=results, alphas=alphas)
misc.toFile(fileName, info)
  Hide output
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-10-1bd6d05d02af> in <module>()
      1 alphas = [-1., -.5, 0., 0.5, 1., 1.5, 2.]
----> 2 fileName = '../data/discriminating_' + info['observer'] + '.pickle'
      3 info['alphas'] = alphas
      4 info['results'] = results
      5 #numpy.savez(fileName, results=results, alphas=alphas)

NameError: name 'info' is not defined

Analysis - version 2

In the second version, we stored more data:

In [11]:
import glob
from psychopy import misc
fig = plt.figure()
ax = fig.add_subplot(111)
for fn in glob.glob('../data/discriminating_v2_*pickle'):
    data = misc.fromFile(fn)
    print fn, results.shape
    print('analyzing results')
    alphas = np.array(data['alphas'])
    print ' alphas = ', alphas
    print '# of trials :', np.abs(data['results']).sum(axis=1)
    print 'average results: ', (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5)
    width = (alphas.max()-alphas.min())/len(alphas)
    ax.bar(data['alphas'] - width/2, (data['results'].sum(axis=1)/np.abs(data['results']).sum(axis=1)*.5+.5), width=width, alpha=.3 )
#    for i_trial in xrange(data['results'].shape[1]):
#        ax.scatter(data['alphas'][data['results']
ax.set_xlabel('alpha')
ax.set_ylabel('% correct')
  Hide output
../data/discriminating_v2_lup_Feb_07_1409.pickle (7, 50)
analyzing results
 alphas =  [-1.  -0.5  0.   0.5  1.   1.5  2. ]
# of trials : [ 10.   4.   5.  11.   6.  11.   3.]
average results:  [ 0.7    1.     0.4    0.545  0.833  0.273  1.   ]
../data/discriminating_v2_lup_Feb_07_1434.pickle (7, 50)
analyzing results
 alphas =  [-1.  -0.5  0.   0.5  1.   1.5  2. ]
# of trials : [  4.   3.  10.  12.   8.   8.   5.]
average results:  [ 0.75   0.667  0.8    0.5    0.125  0.25   0.4  ]

Out[11]:
<matplotlib.text.Text at 0x10d17a610>
In [12]:
# TODO : correct for the number of trials / make global reading / make a new draw of alphas
# TODO : test f_0
  Hide output
In []:
  Hide output

Testing basic functions of Motion Clouds

Motion Clouds: raw principles

Motion Clouds are synthesized textures which aim at having similar characteristics as natural images but with controlled parameters. There are many ways to achieve these results and this notebook aims at showing that different procedures from different communities (neurioscience, modelling, computer vision, ...) may produce similar results.

In [1]:
import numpy as np
np.set_printoptions(precision=3, suppress=True)
import pylab
import matplotlib.pyplot as plt
%matplotlib inline
  Hide output

Using Fourier ("official Motion Clouds")

In [2]:
import sys
sys.path.append('..')
import MotionClouds as mc
fx, fy, ft = mc.get_grids(mc.N_X, mc.N_Y, mc.N_frame)
  Hide output

Using mixtures of images

In [15]:
import scipy.misc
lena = scipy.misc.lena() * 1.
lena -= lena.mean()
lena /= lena.std()
print lena.shape
  Hide output
(512, 512)

In [16]:
plt.imshow(lena, cmap=plt.cm.gray)
  Hide output
Out[16]:
<matplotlib.image.AxesImage at 0x19967c3d0>
In [17]:
lena.shape
  Hide output
Out[17]:
(512, 512)
In [19]:
lena[0, :]
  Hide output
Out[19]:
array([ 0.793,  0.793,  0.793,  0.772,  0.793,  0.689,  0.814,  0.772,
        0.856,  0.772,  0.793,  0.751,  0.647,  0.814,  0.751,  0.647,
        0.689,  0.668,  0.772,  0.772,  0.626,  0.668,  0.626,  0.689,
        0.647,  0.689,  0.647,  0.584,  0.668,  0.626,  0.626,  0.668,
        0.626,  0.709,  0.647,  0.751,  0.709,  0.898,  0.751,  0.877,
        0.877,  0.856,  0.877,  1.002,  0.981,  1.065,  1.023,  0.96 ,
        1.002,  1.002,  0.898,  1.065,  0.918,  0.898,  0.793,  0.772,
        0.501,  0.501,  0.647,  0.354,  0.124, -0.105, -0.126, -0.377,
       -0.565, -0.565, -0.628, -0.67 , -0.753, -0.565, -0.461, -0.586,
       -0.44 , -0.502, -0.419, -0.398, -0.398, -0.398, -0.419, -0.294,
       -0.294, -0.335, -0.356, -0.398, -0.419, -0.314, -0.314, -0.314,
       -0.335, -0.356, -0.356, -0.335, -0.294, -0.314, -0.356, -0.314,
       -0.314, -0.294, -0.294, -0.335, -0.419, -0.377, -0.377, -0.335,
       -0.314, -0.273, -0.147, -0.294, -0.231, -0.126, -0.064, -0.105,
       -0.252, -0.043, -0.064, -0.043,  0.02 , -0.043, -0.064, -0.001,
       -0.001,  0.124,  0.104, -0.022, -0.043,  0.062,  0.166,  0.145,
        0.166,  0.104,  0.229,  0.104,  0.083,  0.145,  0.145,  0.062,
        0.229,  0.25 ,  0.166,  0.208,  0.145,  0.145,  0.083,  0.124,
        0.208,  0.041,  0.104,  0.041,  0.187,  0.124,  0.104,  0.145,
        0.083,  0.166,  0.104,  0.083,  0.229,  0.187,  0.25 ,  0.208,
        0.083,  0.208,  0.124,  0.229,  0.187,  0.187,  0.208,  0.229,
        0.166,  0.229,  0.145,  0.25 ,  0.208,  0.166,  0.187,  0.187,
        0.25 ,  0.187,  0.208,  0.208,  0.166,  0.187,  0.145,  0.166,
        0.187,  0.229,  0.271,  0.229,  0.187,  0.145,  0.187,  0.145,
        0.229,  0.25 ,  0.166,  0.229,  0.25 ,  0.271,  0.187,  0.292,
        0.25 ,  0.208,  0.187,  0.229,  0.229,  0.25 ,  0.166,  0.187,
        0.333,  0.166,  0.145,  0.312,  0.208,  0.145,  0.104,  0.166,
        0.166,  0.062,  0.145,  0.166,  0.166,  0.145,  0.104,  0.166,
        0.145,  0.083,  0.229,  0.187,  0.187,  0.229,  0.145,  0.166,
        0.166,  0.166,  0.145,  0.187,  0.124,  0.145,  0.145,  0.187,
        0.25 ,  0.208,  0.166,  0.229,  0.166,  0.229,  0.208,  0.292,
        0.166,  0.166,  0.187,  0.208,  0.104,  0.104,  0.187,  0.292,
        0.292,  0.375,  0.208,  0.166,  0.166,  0.104,  0.083,  0.083,
        0.166,  0.208,  0.104,  0.145,  0.208,  0.187,  0.145,  0.083,
        0.104,  0.062,  0.145,  0.166,  0.041,  0.124,  0.229,  0.208,
        0.145,  0.062,  0.041,  0.145,  0.062,  0.145,  0.208,  0.062,
        0.041,  0.166,  0.145,  0.083,  0.145,  0.062,  0.145,  0.02 ,
        0.062,  0.124,  0.02 ,  0.062,  0.062,  0.166,  0.041,  0.104,
       -0.064, -0.001, -0.001, -0.064, -0.105, -0.043, -0.064, -0.105,
       -0.147, -0.189, -0.252, -0.231, -0.189, -0.398, -0.461, -0.44 ,
       -0.419, -0.252, -0.168, -0.126, -0.043,  0.166,  0.229,  0.312,
        0.333,  0.417,  0.459,  0.542,  0.605,  0.584,  0.605,  0.73 ,
        0.73 ,  0.73 ,  0.814,  0.709,  0.793,  0.772,  0.647,  0.647,
        0.501,  0.542,  0.605,  0.542,  0.626,  0.563,  0.563,  0.647,
        0.563,  0.689,  0.647,  0.584,  0.689,  0.668,  0.626,  0.626,
        0.73 ,  0.709,  0.647,  0.689,  0.626,  0.605,  0.563,  0.689,
        0.605,  0.563,  0.626,  0.584,  0.626,  0.605,  0.584,  0.605,
        0.605,  0.563,  0.563,  0.709,  0.626,  0.626,  0.584,  0.584,
        0.605,  0.689,  0.626,  0.709,  0.709,  0.689,  0.626,  0.626,
        0.668,  0.73 ,  0.73 ,  0.73 ,  0.689,  0.668,  0.668,  0.584,
        0.647,  0.647,  0.709,  1.107,  1.42 ,  1.545,  1.713,  1.754,
        1.859,  1.838,  1.859,  1.921,  1.963,  1.963,  1.963,  1.963,
        2.005,  1.921,  1.817,  1.671,  1.42 ,  1.002,  0.584, -0.043,
       -0.377, -0.398, -0.482, -0.44 , -0.419, -0.273, -0.398, -0.231,
       -0.168, -0.147, -0.105, -0.147, -0.105, -0.126, -0.064, -0.022,
        0.062, -0.064, -0.105, -0.147, -0.085,  0.02 , -0.064, -0.126,
       -0.022, -0.043,  0.02 ,  0.083, -0.064, -0.105, -0.022, -0.126,
       -0.085, -0.022, -0.022, -0.105, -0.043, -0.064, -0.064, -0.022,
       -0.043, -0.001,  0.104, -0.126, -0.001,  0.041, -0.001,  0.041,
        0.041,  0.062, -0.022, -0.001, -0.001,  0.02 , -0.126, -0.126,
        0.083, -0.043,  0.041,  0.041,  0.041,  0.124, -0.001,  0.041,
        0.02 , -0.064,  0.104, -0.001,  0.041,  0.062, -0.105, -0.022,
       -0.21 , -0.064, -0.168, -0.105, -0.126, -0.189, -0.085, -0.001,
        0.292,  0.835,  0.918,  0.96 ,  0.981,  0.96 ,  0.647,  0.083])
In [5]:
def noise(image=lena):
    for axis in [0, 1]:
        image = np.roll(image, np.random.randint(image.shape[axis]), axis=axis)
    return image
  Hide output
In [6]:
plt.imshow(noise(), cmap=plt.cm.gray)
  Hide output
Out[6]:
<matplotlib.image.AxesImage at 0x117793790>
In [7]:
plt.imshow(noise(), cmap=plt.cm.gray)
  Hide output
Out[7]:
<matplotlib.image.AxesImage at 0x117c95890>

Using ARMA processes

Now, we define the ARMA process as an averaging process with a certain time constant $\tau=30.$ (in frames).

In [8]:
def ARMA(image, tau=30.):
    image = (1 - 1/tau)* image + 1/tau * noise()
    return image
  Hide output

initializing

In [9]:
image = ARMA(lena)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[9]:
<matplotlib.image.AxesImage at 0x118193b50>
In [10]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[10]:
<matplotlib.image.AxesImage at 0x11820bcd0>
In [11]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[11]:
<matplotlib.image.AxesImage at 0x118bccd10>
In [12]:
for _ in range(1000): image = ARMA(image)
plt.imshow(image, cmap=plt.cm.gray)
  Hide output
Out[12]:
<matplotlib.image.AxesImage at 0x11910b8d0>
In [13]:
%cd..
N_frame = 1024
z = np.zeros((lena.shape[0], lena.shape[1], N_frame))
z[:, :, 0] = image
for i_frame in range(1, N_frame): 
    z[:, :, i_frame] = ARMA(z[:, :, i_frame-1])
mc.anim_save(.5 + .5*z, filename='results/arma')
  Hide output
[  0 %] elapsed[sec]: 0.000
/Users/lolo/Dropbox/science/MotionClouds
Saving sequence results/arma.webm
[51150 %] elapsed[sec]: 57.826 | ETA[sec]: -57.713 

Title: 
                      Started: 09/23/2014 11:58:20
                      Finished: 09/23/2014 11:59:18
                      Total time elapsed: 0.000 sec

In [14]:
mc.notebook = True
mc.in_show_video(name='arma')
  Hide output
arma

    Share